The aim of these exercises is to continue making sure you’re comfortable with handling multivariate data. In this chapter’s, you’ll focus on analyses based on dissimilarities/distances, including fitting linear models to these kinds of response variables.
For each of the examples below, you should follow the sequence we’ve used previously, as far as it’s sensible:
What is the biological question?
Is the predictor continuous or categorical?
Write out the linear model corresponding to this question.
What distribution do you expect the response variable to follow?
What are the assumptions behind the statistical model you’ll fit?
Fit the model
How will you assess whether the model fits well?
Can you detect an effect of the predictors?
How do you measure the effect?
What do you conclude (including any cautions)
Dixon et al. (2018) focused on assemblages of reptiles in forests and woodlands of southeastern Australia, with a particular interest in relationships between these assemblages and the fire history of the landscape. They identified 81 sites that varied in time since fire, from 6 months to >96 y. Rather than a continuum, three categories of time since fire were used, 0.5-2y, 6-12y, and >96y. Sites were also classified according to habitat.
Reptile assemblages were sampled with a range of methods, including visual surveys and camera traps, to give counts of 20 reptiles, whose abundance ranged from 0-1 to 0-126, depending on species.
Data are available from dryad, as Rep_abund.csv. You’ll want to focus on the columns with reptile numbers, and the one with the fire history (tsf). For convenience, we’ve extracted those data as dixonbiota.csv
Start by having a quick look at the data.
df <- read_csv("data/dixonbiota.csv")
Rows: 79 Columns: 22── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): site, tsf
dbl (20): adup, aplat, amur, amac, aram, dcor, ecun, esax, eul, htal, ldel, lgui, lwhi, ppor, pent, pspe, ptex, rdie, tnig, vros
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(df, 10)
The file is pretty straightforward, with tsf being the fire category, and columns 3-22 each representing one reptile taxon.
Repeat the preceding analysis after creating a new data file that is presence-absence only.
df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
#Run 3d first
df1s.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.1404128
Run 1 stress 0.1480196
Run 2 stress 0.1405858
... Procrustes: rmse 0.01565156 max resid 0.08607928
Run 3 stress 0.1485349
Run 4 stress 0.1417498
Run 5 stress 0.1404104
... New best solution
... Procrustes: rmse 0.0023304 max resid 0.01808758
Run 6 stress 0.1404131
... Procrustes: rmse 0.002334262 max resid 0.01834106
Run 7 stress 0.1415289
Run 8 stress 0.140411
... Procrustes: rmse 0.0007835425 max resid 0.005200482
... Similar to previous best
Run 9 stress 0.1486406
Run 10 stress 0.1404103
... New best solution
... Procrustes: rmse 0.001207819 max resid 0.009469643
... Similar to previous best
Run 11 stress 0.1404117
... Procrustes: rmse 0.000916645 max resid 0.006636852
... Similar to previous best
Run 12 stress 0.1405581
... Procrustes: rmse 0.01501831 max resid 0.08924602
Run 13 stress 0.143846
Run 14 stress 0.1495992
Run 15 stress 0.140558
... Procrustes: rmse 0.01458579 max resid 0.08787773
Run 16 stress 0.1404148
... Procrustes: rmse 0.0015373 max resid 0.01187026
Run 17 stress 0.1448137
Run 18 stress 0.1443774
Run 19 stress 0.140411
... Procrustes: rmse 0.001777414 max resid 0.01404256
Run 20 stress 0.1500017
Run 21 stress 0.1485371
Run 22 stress 0.1454737
Run 23 stress 0.140411
... Procrustes: rmse 0.00189255 max resid 0.01470197
Run 24 stress 0.1496201
Run 25 stress 0.1405615
... Procrustes: rmse 0.01445625 max resid 0.08690181
Run 26 stress 0.1477962
Run 27 stress 0.1478653
Run 28 stress 0.1488089
Run 29 stress 0.1444777
Run 30 stress 0.1405647
... Procrustes: rmse 0.01491986 max resid 0.08994996
Run 31 stress 0.1405594
... Procrustes: rmse 0.01486184 max resid 0.08799154
Run 32 stress 0.1405612
... Procrustes: rmse 0.01487556 max resid 0.0876917
Run 33 stress 0.140559
... Procrustes: rmse 0.01485206 max resid 0.0880716
Run 34 stress 0.1425376
Run 35 stress 0.140417
... Procrustes: rmse 0.001866647 max resid 0.01406319
Run 36 stress 0.1418405
Run 37 stress 0.1422689
Run 38 stress 0.1486242
Run 39 stress 0.1413517
Run 40 stress 0.1468941
*** Best solution repeated 2 times
stressplot(df1s.mds3, main="Shepard plot")
df1s.mds3
Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.jac
Distance: binary jaccard
Dimensions: 3
Stress: 0.1404103
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 10 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
# Now look at 2d
df1s.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.2009509
Run 1 stress 0.2079375
Run 2 stress 0.2224574
Run 3 stress 0.2162073
Run 4 stress 0.225869
Run 5 stress 0.2284785
Run 6 stress 0.2305984
Run 7 stress 0.2199417
Run 8 stress 0.202096
Run 9 stress 0.2011626
... Procrustes: rmse 0.006551419 max resid 0.03916879
Run 10 stress 0.2104293
Run 11 stress 0.2015109
Run 12 stress 0.2009556
... Procrustes: rmse 0.001703517 max resid 0.01335007
Run 13 stress 0.2135392
Run 14 stress 0.2246277
Run 15 stress 0.2240908
Run 16 stress 0.2250626
Run 17 stress 0.200954
... Procrustes: rmse 0.002998174 max resid 0.01914574
Run 18 stress 0.2095741
Run 19 stress 0.2186864
Run 20 stress 0.203307
Run 21 stress 0.2081138
Run 22 stress 0.2212435
Run 23 stress 0.2132584
Run 24 stress 0.2055683
Run 25 stress 0.2019164
Run 26 stress 0.2108398
Run 27 stress 0.2151842
Run 28 stress 0.2009453
... New best solution
... Procrustes: rmse 0.002748903 max resid 0.0197823
Run 29 stress 0.206872
Run 30 stress 0.2224539
Run 31 stress 0.2031111
Run 32 stress 0.2064081
Run 33 stress 0.2149136
Run 34 stress 0.2009451
... New best solution
... Procrustes: rmse 0.002076934 max resid 0.0168433
Run 35 stress 0.2111817
Run 36 stress 0.2164562
Run 37 stress 0.2027405
Run 38 stress 0.2161414
Run 39 stress 0.2242858
Run 40 stress 0.2071548
Run 41 stress 0.203393
Run 42 stress 0.2266633
Run 43 stress 0.2121352
Run 44 stress 0.2175283
Run 45 stress 0.2104325
Run 46 stress 0.2029998
Run 47 stress 0.2159185
Run 48 stress 0.2042398
Run 49 stress 0.2134569
Run 50 stress 0.2145411
Run 51 stress 0.2081612
Run 52 stress 0.223412
Run 53 stress 0.2327611
Run 54 stress 0.2199861
Run 55 stress 0.2115448
Run 56 stress 0.209779
Run 57 stress 0.2019152
Run 58 stress 0.2236233
Run 59 stress 0.2015106
Run 60 stress 0.2036071
Run 61 stress 0.2118656
Run 62 stress 0.2120415
Run 63 stress 0.2222671
Run 64 stress 0.2176459
Run 65 stress 0.2133805
Run 66 stress 0.2009556
... Procrustes: rmse 0.003462696 max resid 0.0201114
Run 67 stress 0.201907
Run 68 stress 0.2009456
... Procrustes: rmse 0.0003121967 max resid 0.001856644
... Similar to previous best
*** Best solution repeated 1 times
stressplot(df1s.mds2, main="Shepard plot")
df1s.mds2
Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.jac
Distance: binary jaccard
Dimensions: 2
Stress: 0.2009451
Stress type 1, weak ties
Best solution was repeated 1 time in 68 tries
The best solution was from try 34 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
Marginal decision here - stress just over 0.2 for k=2, so might be OK. Stress good for k=3, but 2-d MDS usually better for reporting or showing to audience
a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:2)],a) #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=tsf, ) )+
geom_point()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=sym3,
name="TSF",
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = tsf, ) )+
geom_point()+
labs(y="MDS3", x="MDS1")+
scale_shape_manual(values=sym3,
name="TSF",
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()
Unburned sites separate very clearly on MDS1. Pattern clearer than with abundances included
Do permanova including tsf
df1.ado <- adonis2(df1.jac~tsf,data=df,permutations=999)
print(df1.ado)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = df1.jac ~ tsf, data = df, permutations = 999)
Df SumOfSqs R2 F Pr(>F)
tsf 2 3.7853 0.20678 9.9057 0.001 ***
Residual 76 14.5210 0.79322
Total 78 18.3062 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Clear fire signal
Dixon and her colleagues also recorded a range of habitat variables, which we’ve collected for you in dixonenv.csv
dixonenv <- read_csv("data/dixonenv.csv")
Rows: 79 Columns: 14── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (4): site, tsf, veg, aspect
dbl (10): lcwd, shrcov, grcov, litcov, litdep, rocks, elev, warm, cold, twi
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(dixonenv, 10)
They recorded Coarse Woody Debris, which was long-transformed, litter cover (litcov), % cover of groundcover (grcov), and % cover of shrubs (shrcov), and rock cover. In the original paper, Dixon et al. were comfortable that these predictors weren’t correlated.
dixonbiotamv <- mvabund(dixonbiota[,-(1:2)])
dixonbiotamv.mv <- manyglm(dixonbiotamv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="poisson")
plot(dixonbiotamv.mv)
# still hetero vars so use -ve binomial
dixonbiotamv1.mv <- manyglm(dixonbiotamv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="negative_binomial")
plot(dixonbiotamv1.mv)
anova(dixonbiotamv1.mv)
Time elapsed: 0 hr 1 min 39 sec
Analysis of Deviance Table
Model: dixonbiotamv ~ dixonenv$tsf + dixonenv$lcwd + dixonenv$litcov + dixonenv$grcov + dixonenv$shrcov + dixonenv$lrocks
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 78
dixonenv$tsf 76 2 331.6 0.001 ***
dixonenv$lcwd 75 1 61.5 0.001 ***
dixonenv$litcov 74 1 25.9 0.265
dixonenv$grcov 73 1 36.4 0.046 *
dixonenv$shrcov 72 1 21.5 0.607
dixonenv$lrocks 71 1 58.0 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 iterations via PIT-trap resampling.
Let’s return to the Hutto and Barrett (2021) example from the Chapter 15 exercises. There, we used similarity-based analyses to explore the relationship between frog assemblages in ponds and the degree of urbanization surrounding. This seems a question that could just as easily be examined using distances or dissimilarities.
Return to the data and assess whether frog assemblages (using the standardized abundance scale or presence-absence) differ with urbanization.
df <- read_csv("data/huttoamph.csv")
Rows: 50 Columns: 14── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): type
dbl (13): Site, ACRE, BAME, BFOW, GCAR, HCIN, HVER, LCAT, LCLA, LPAL, LSPH, PCRU, PFER
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(df,10)
#ponds 9 and 40 had no frogs. Remove for distance analysis
df <- df %>%
filter(Site!= 9 & Site !=40)
df
The explanation of the variables is in the previous chapter’s exercises.
df1 <- df[,-(1:2)] #Stripped back file without labels, leaving only
df1.bc <- vegdist(df1,'bray')
#Run 2d first
df1.mds2 <- metaMDS(df1.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.196708
Run 1 stress 0.1877626
... New best solution
... Procrustes: rmse 0.06205021 max resid 0.2289634
Run 2 stress 0.1888251
Run 3 stress 0.1920872
Run 4 stress 0.1853642
... New best solution
... Procrustes: rmse 0.07936891 max resid 0.4999309
Run 5 stress 0.1855846
... Procrustes: rmse 0.01953946 max resid 0.1056279
Run 6 stress 0.1984394
Run 7 stress 0.1862238
Run 8 stress 0.200767
Run 9 stress 0.1961235
Run 10 stress 0.1853641
... New best solution
... Procrustes: rmse 9.823368e-05 max resid 0.0004897303
... Similar to previous best
Run 11 stress 0.204919
Run 12 stress 0.2001183
Run 13 stress 0.1875595
Run 14 stress 0.2031338
Run 15 stress 0.1875596
Run 16 stress 0.2001202
Run 17 stress 0.1877626
Run 18 stress 0.1968146
Run 19 stress 0.2191092
Run 20 stress 0.1862237
Run 21 stress 0.1862237
Run 22 stress 0.2016092
Run 23 stress 0.1919379
Run 24 stress 0.1921465
Run 25 stress 0.1875596
Run 26 stress 0.1862241
Run 27 stress 0.1888273
Run 28 stress 0.1958875
Run 29 stress 0.2103827
Run 30 stress 0.196403
Run 31 stress 0.1920872
Run 32 stress 0.1853644
... Procrustes: rmse 0.0001814753 max resid 0.0007769668
... Similar to previous best
Run 33 stress 0.1929083
Run 34 stress 0.1948815
Run 35 stress 0.2023845
Run 36 stress 0.1961143
Run 37 stress 0.1978435
Run 38 stress 0.214542
Run 39 stress 0.1970125
Run 40 stress 0.2055679
*** Best solution repeated 2 times
stressplot(df1.mds2, main="Shepard plot")
df1.mds2
Call:
metaMDS(comm = df1.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.bc
Distance: bray
Dimensions: 2
Stress: 0.1853641
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 10 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
# Now look at 3d
df.mds3 <- metaMDS(df1.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.129187
Run 1 stress 0.1388768
Run 2 stress 0.1382418
Run 3 stress 0.1383885
Run 4 stress 0.129186
... New best solution
... Procrustes: rmse 0.001355541 max resid 0.007232405
... Similar to previous best
Run 5 stress 0.129186
... Procrustes: rmse 0.0005937864 max resid 0.001870774
... Similar to previous best
Run 6 stress 0.1291858
... New best solution
... Procrustes: rmse 0.0008376747 max resid 0.004337427
... Similar to previous best
Run 7 stress 0.129186
... Procrustes: rmse 0.0007498725 max resid 0.00420023
... Similar to previous best
Run 8 stress 0.1291864
... Procrustes: rmse 0.0002763668 max resid 0.0009518002
... Similar to previous best
Run 9 stress 0.1383896
Run 10 stress 0.1384317
Run 11 stress 0.1289614
... New best solution
... Procrustes: rmse 0.0109308 max resid 0.05605748
Run 12 stress 0.128962
... Procrustes: rmse 0.0009039712 max resid 0.005251417
... Similar to previous best
Run 13 stress 0.1383929
Run 14 stress 0.1402552
Run 15 stress 0.128962
... Procrustes: rmse 0.0009246251 max resid 0.005329771
... Similar to previous best
Run 16 stress 0.1291859
... Procrustes: rmse 0.01058901 max resid 0.0547769
Run 17 stress 0.1289615
... Procrustes: rmse 0.000419003 max resid 0.002034791
... Similar to previous best
Run 18 stress 0.1380889
Run 19 stress 0.129186
... Procrustes: rmse 0.01092673 max resid 0.05630155
Run 20 stress 0.1383886
Run 21 stress 0.1289617
... Procrustes: rmse 0.0002465603 max resid 0.0007331739
... Similar to previous best
Run 22 stress 0.1291604
... Procrustes: rmse 0.009890855 max resid 0.05184425
Run 23 stress 0.1291601
... Procrustes: rmse 0.009490803 max resid 0.05078321
Run 24 stress 0.1291855
... Procrustes: rmse 0.01066625 max resid 0.05496251
Run 25 stress 0.1391225
Run 26 stress 0.129186
... Procrustes: rmse 0.01099671 max resid 0.05629543
Run 27 stress 0.1291866
... Procrustes: rmse 0.01114331 max resid 0.05711187
Run 28 stress 0.1384282
Run 29 stress 0.1291861
... Procrustes: rmse 0.01101164 max resid 0.0563452
Run 30 stress 0.1383979
Run 31 stress 0.1289613
... New best solution
... Procrustes: rmse 0.0004223583 max resid 0.001654659
... Similar to previous best
Run 32 stress 0.1380874
Run 33 stress 0.1291858
... Procrustes: rmse 0.01077276 max resid 0.05558138
Run 34 stress 0.1291859
... Procrustes: rmse 0.01081878 max resid 0.05582879
Run 35 stress 0.1291858
... Procrustes: rmse 0.01042548 max resid 0.05433355
Run 36 stress 0.1385874
Run 37 stress 0.1291863
... Procrustes: rmse 0.01092361 max resid 0.05638233
Run 38 stress 0.1308828
Run 39 stress 0.1384275
Run 40 stress 0.1388191
*** Best solution repeated 1 times
stressplot(df.mds3, main="Shepard plot")
df.mds3
Call:
metaMDS(comm = df1.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.bc
Distance: bray
Dimensions: 3
Stress: 0.1289613
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 31 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
2-d is acceptable
a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a) #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
geom_point()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=sym3,
name="Urbanization",
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)+
theme_qk() + scale_color_uchicago()
df1.ado <- adonis2(df1.bc~type,data=df,permutations=999)
print(df1.ado)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = df1.bc ~ type, data = df, permutations = 999)
Df SumOfSqs R2 F Pr(>F)
type 2 0.912 0.07281 1.7669 0.045 *
Residual 45 11.614 0.92719
Total 47 12.526 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permanova detects a difference; MDS plot suggests low urbanization may be different from the other two, though worth looking at dispersion as well, as spread of this type is less than the other two on plot
Repeat analysis with presence-absence
df1.jac <- vegdist(df1,binary=TRUE,method=‘jaccard’)
#Run 2d first
df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
df1.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.1551612
Run 1 stress 0.1630921
Run 2 stress 0.1622263
Run 3 stress 0.1724054
Run 4 stress 0.1585799
Run 5 stress 0.1761433
Run 6 stress 0.1608258
Run 7 stress 0.1734443
Run 8 stress 0.1621126
Run 9 stress 0.1536788
... New best solution
... Procrustes: rmse 0.01901219 max resid 0.06506965
Run 10 stress 0.1760283
Run 11 stress 0.1609339
Run 12 stress 0.1565712
Run 13 stress 0.1939366
Run 14 stress 0.1843833
Run 15 stress 0.162526
Run 16 stress 0.1622263
Run 17 stress 0.1570453
Run 18 stress 0.1708675
Run 19 stress 0.1537515
... Procrustes: rmse 0.003076067 max resid 0.01474863
Run 20 stress 0.1539129
... Procrustes: rmse 0.00947995 max resid 0.04425278
Run 21 stress 0.1664178
Run 22 stress 0.1795028
Run 23 stress 0.1585801
Run 24 stress 0.1555104
Run 25 stress 0.1708671
Run 26 stress 0.167874
Run 27 stress 0.1673092
Run 28 stress 0.1719174
Run 29 stress 0.1591019
Run 30 stress 0.1711482
Run 31 stress 0.1795847
Run 32 stress 0.1608763
Run 33 stress 0.1721431
Run 34 stress 0.1585798
Run 35 stress 0.1850134
Run 36 stress 0.1607836
Run 37 stress 0.1586731
Run 38 stress 0.1738412
Run 39 stress 0.2004114
Run 40 stress 0.1672269
Run 41 stress 0.1585799
Run 42 stress 0.1591019
Run 43 stress 0.162564
Run 44 stress 0.1620556
Run 45 stress 0.1679859
Run 46 stress 0.171033
Run 47 stress 0.1617404
Run 48 stress 0.1538246
... Procrustes: rmse 0.01103138 max resid 0.05221548
Run 49 stress 0.172142
Run 50 stress 0.1960733
Run 51 stress 0.1731862
Run 52 stress 0.1627813
Run 53 stress 0.1536788
... Procrustes: rmse 0.0001144048 max resid 0.000322326
... Similar to previous best
*** Best solution repeated 1 times
stressplot(df1.mds2, main="Shepard plot")
df1.mds2
Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.jac
Distance: binary jaccard
Dimensions: 2
Stress: 0.1536788
Stress type 1, weak ties
Best solution was repeated 1 time in 53 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
# Now look at 3d
df.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.09466813
Run 1 stress 0.09466802
... New best solution
... Procrustes: rmse 0.0004764601 max resid 0.002102695
... Similar to previous best
Run 2 stress 0.09466791
... New best solution
... Procrustes: rmse 0.0001518385 max resid 0.0004555261
... Similar to previous best
Run 3 stress 0.09466803
... Procrustes: rmse 0.0001108676 max resid 0.0004745956
... Similar to previous best
Run 4 stress 0.1005052
Run 5 stress 0.09466787
... New best solution
... Procrustes: rmse 5.340853e-05 max resid 0.0001992126
... Similar to previous best
Run 6 stress 0.1005044
Run 7 stress 0.1005056
Run 8 stress 0.09466803
... Procrustes: rmse 0.0001288817 max resid 0.0005023949
... Similar to previous best
Run 9 stress 0.1005047
Run 10 stress 0.09466793
... Procrustes: rmse 5.310415e-05 max resid 0.0002709636
... Similar to previous best
Run 11 stress 0.09466796
... Procrustes: rmse 0.0001256892 max resid 0.0003435652
... Similar to previous best
Run 12 stress 0.094668
... Procrustes: rmse 9.745181e-05 max resid 0.0004669281
... Similar to previous best
Run 13 stress 0.09693819
Run 14 stress 0.09694632
Run 15 stress 0.1005051
Run 16 stress 0.09466833
... Procrustes: rmse 0.0002133226 max resid 0.001028399
... Similar to previous best
Run 17 stress 0.09466821
... Procrustes: rmse 0.0002054658 max resid 0.0009036094
... Similar to previous best
Run 18 stress 0.09466811
... Procrustes: rmse 0.0001423043 max resid 0.0004843394
... Similar to previous best
Run 19 stress 0.09466805
... Procrustes: rmse 0.0003941737 max resid 0.002017784
... Similar to previous best
Run 20 stress 0.109714
Run 21 stress 0.1005051
Run 22 stress 0.09466816
... Procrustes: rmse 0.000195071 max resid 0.0006552548
... Similar to previous best
Run 23 stress 0.09466783
... New best solution
... Procrustes: rmse 6.655458e-05 max resid 0.0002929959
... Similar to previous best
Run 24 stress 0.09466805
... Procrustes: rmse 0.0001689595 max resid 0.0006247545
... Similar to previous best
Run 25 stress 0.1005048
Run 26 stress 0.09466801
... Procrustes: rmse 0.0001773767 max resid 0.0009243069
... Similar to previous best
Run 27 stress 0.09694616
Run 28 stress 0.1004898
Run 29 stress 0.09466783
... New best solution
... Procrustes: rmse 0.0001519669 max resid 0.0004171994
... Similar to previous best
Run 30 stress 0.1005053
Run 31 stress 0.1097146
Run 32 stress 0.1004899
Run 33 stress 0.09466786
... Procrustes: rmse 0.0001828158 max resid 0.000484661
... Similar to previous best
Run 34 stress 0.1090648
Run 35 stress 0.09466812
... Procrustes: rmse 0.0002666448 max resid 0.0009620816
... Similar to previous best
Run 36 stress 0.09694549
Run 37 stress 0.09694624
Run 38 stress 0.094668
... Procrustes: rmse 0.0002647032 max resid 0.001383475
... Similar to previous best
Run 39 stress 0.09466798
... Procrustes: rmse 0.0002442504 max resid 0.0007524015
... Similar to previous best
Run 40 stress 0.09466802
... Procrustes: rmse 0.000254564 max resid 0.0008115543
... Similar to previous best
*** Best solution repeated 6 times
stressplot(df.mds3, main="Shepard plot")
df.mds3
Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1.jac
Distance: binary jaccard
Dimensions: 3
Stress: 0.09466783
Stress type 1, weak ties
Best solution was repeated 6 times in 40 tries
The best solution was from try 29 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
Again, stress OK with 2D, good with 3D. Look at 2D first
a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a) #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
geom_point()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=sym3,
name="Urbanization",
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)+
theme_qk() + scale_color_uchicago()
df1.ado <- adonis2(df1.jac~type,data=df,permutations=999)
print(df1.ado)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = df1.jac ~ type, data = df, permutations = 999)
Df SumOfSqs R2 F Pr(>F)
type 2 0.4525 0.04389 1.0328 0.403
Residual 45 9.8587 0.95611
Total 47 10.3112 1.00000
No separation seen on Presence-Absence
Griffen et al. (2019) examined the oral microbiome of humans, with a focus on the effects of HIV and AntiRetroviral therapy (ART) on this microbiome. They described the microbiome of 341 patients, who fell into one of two categories: HIV-, and HIV+ with ART. They also matched their samples as far as possible for sex, along with other conditions that can influence oral microbiomes (Candida infection, current smoking). We’ll use their records for these other categories as well. There were more HIV+ than - subjects, but within these groups, approximately equal numbers of two sexes, two Candida infection status, and two smoking categories.
The bacteriome was recorded separately using 16S RNA sequencing, which overed just over 600 taxa.
The data are available from their paper, as two Excel files in the supplementary information. The metadata gives HIV status, and a raft of demographic information. We’ll just work with HIV, sex, Candida, and current smoking. A second sheet has the bacteriome. The code chunk below reads this file (from a Downloads folder - you’ll need to modify the file location). It also screens the file for any bacterial taxa that were not present in this particular study, which reduces the bacterial taxa to 599.
library(readxl)
#Read in metadata. Lots of information we don't need, so a couple of iterations to get down to a few columns we want - HIV status, candida, current smoking, gender
df <- read_excel("~/Downloads/41598_2019_55703_MOESM3_ESM.xlsx",
range = "A1:S342", na = "NA")
df<-select(df, HIV, candida, gender, smoking_current)
head(df)
#df1 is the bacterial data
df1 <- read_excel("~/Downloads/41598_2019_55703_MOESM2_ESM.xlsx")
New names:
head(df1,10)
df1 <- df1 %>%
select(where( ~ is.numeric(.x) && sum(.x) != 0)) #Drops any bacterial taxon not recorded in this study (i.e. where it's all zeroes)
df1 <- df1[,-1] #remove col1, which is sample ID
Do these factors act independently on the bacteriome?
Which effects are largest (use some graphical methods)?
First steps: Think about standardization and which distance measure you’ll use. Bray-Curtis seems common for these kind of data, and and counts of different taxa vary widely.
df1s <- wisconsin(df1) #Common to do some standardisation, and counts of different taxa vary widely
df1s.bc <- vegdist(df1s,'bray') #With standardisation (wisconsin)
Run factorial model using permanova, based on B-C
df1.ado <- adonis2(df1s.bc~HIV*candida*gender*smoking_current,data=df,permutations=999)
print(df1.ado)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = df1s.bc ~ HIV * candida * gender * smoking_current, data = df, permutations = 999)
Df SumOfSqs R2 F Pr(>F)
HIV 1 1.369 0.01250 4.3776 0.001 ***
candida 1 1.461 0.01333 4.6694 0.001 ***
gender 1 0.492 0.00449 1.5737 0.009 **
smoking_current 1 1.111 0.01014 3.5513 0.001 ***
HIV:candida 1 0.318 0.00290 1.0151 0.401
HIV:gender 1 0.354 0.00324 1.1332 0.204
candida:gender 1 0.257 0.00234 0.8202 0.861
HIV:smoking_current 1 0.269 0.00245 0.8584 0.765
candida:smoking_current 1 0.344 0.00314 1.0989 0.243
gender:smoking_current 1 0.394 0.00360 1.2602 0.098 .
HIV:candida:gender 1 0.349 0.00319 1.1164 0.203
HIV:candida:smoking_current 1 0.280 0.00255 0.8943 0.707
HIV:gender:smoking_current 1 0.265 0.00242 0.8463 0.793
candida:gender:smoking_current 1 0.254 0.00232 0.8110 0.841
HIV:candida:gender:smoking_current 1 0.362 0.00330 1.1562 0.213
Residual 325 101.664 0.92809
Total 340 109.541 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Model shows main effects, but no interactions important (or at least, detected)
Could you fit a simpler model to the data?
Run simpler model
df2.ado <- adonis2(df1s.bc~HIV+candida+gender+smoking_current,data=df,permutations=999)
print(df2.ado)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = df1s.bc ~ HIV + candida + gender + smoking_current, data = df, permutations = 999)
Df SumOfSqs R2 F Pr(>F)
HIV 1 1.369 0.01250 4.3775 0.001 ***
candida 1 1.461 0.01333 4.6693 0.001 ***
gender 1 0.492 0.00449 1.5736 0.013 *
smoking_current 1 1.111 0.01014 3.5512 0.001 ***
Residual 336 105.108 0.95953
Total 340 109.541 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
No change; not too surprising, as all four main effects were detected easily.
Now go to data visualization
#Run 3d first
df1s.mds3 <- metaMDS(df1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.1573956
Run 1 stress 0.1579391
Run 2 stress 0.1572758
... New best solution
... Procrustes: rmse 0.01025549 max resid 0.1525195
Run 3 stress 0.1579179
Run 4 stress 0.157493
... Procrustes: rmse 0.007283638 max resid 0.1009598
Run 5 stress 0.1578733
Run 6 stress 0.1593345
Run 7 stress 0.1596534
Run 8 stress 0.1578064
Run 9 stress 0.1576408
... Procrustes: rmse 0.005563033 max resid 0.08072257
Run 10 stress 0.1574589
... Procrustes: rmse 0.007053274 max resid 0.1007079
Run 11 stress 0.1574113
... Procrustes: rmse 0.008390411 max resid 0.1530435
Run 12 stress 0.1584759
Run 13 stress 0.1590917
Run 14 stress 0.1579313
Run 15 stress 0.1577686
... Procrustes: rmse 0.00731071 max resid 0.08902012
Run 16 stress 0.1579681
Run 17 stress 0.157889
Run 18 stress 0.1576122
... Procrustes: rmse 0.01046473 max resid 0.1079021
Run 19 stress 0.1575193
... Procrustes: rmse 0.01344117 max resid 0.1517403
Run 20 stress 0.1590183
Run 21 stress 0.1572873
... Procrustes: rmse 0.005938186 max resid 0.1084724
Run 22 stress 0.1579056
Run 23 stress 0.1578609
Run 24 stress 0.1575202
... Procrustes: rmse 0.01348549 max resid 0.1517058
Run 25 stress 0.1580688
Run 26 stress 0.1573537
... Procrustes: rmse 0.005408552 max resid 0.09851134
Run 27 stress 0.1584273
Run 28 stress 0.1575528
... Procrustes: rmse 0.007336861 max resid 0.1083012
Run 29 stress 0.1591385
Run 30 stress 0.1575529
... Procrustes: rmse 0.007325558 max resid 0.1082856
Run 31 stress 0.1584616
Run 32 stress 0.157776
Run 33 stress 0.1579148
Run 34 stress 0.1590267
Run 35 stress 0.1574802
... Procrustes: rmse 0.01063831 max resid 0.1527981
Run 36 stress 0.1574825
... Procrustes: rmse 0.0125045 max resid 0.1520643
Run 37 stress 0.158388
Run 38 stress 0.1581063
Run 39 stress 0.1586636
Run 40 stress 0.159344
Run 41 stress 0.1575585
... Procrustes: rmse 0.007481761 max resid 0.1081395
Run 42 stress 0.1580441
Run 43 stress 0.1583809
Run 44 stress 0.1586431
Run 45 stress 0.1575179
... Procrustes: rmse 0.009346751 max resid 0.1082067
Run 46 stress 0.1585796
Run 47 stress 0.1578686
Run 48 stress 0.1595925
Run 49 stress 0.1579636
Run 50 stress 0.1572746
... New best solution
... Procrustes: rmse 0.0005494158 max resid 0.00824335
... Similar to previous best
*** Best solution repeated 1 times
stressplot(df1s.mds3, main="Shepard plot")
df1s.mds3
Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1s.bc
Distance: bray
Dimensions: 3
Stress: 0.1572746
Stress type 1, weak ties
Best solution was repeated 1 time in 50 tries
The best solution was from try 50 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
# Now look at 2d
df1s.mds2 <- metaMDS(df1s.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)
Run 0 stress 0.2249552
Run 1 stress 0.2254408
... Procrustes: rmse 0.01613467 max resid 0.2308209
Run 2 stress 0.2253412
... Procrustes: rmse 0.01606138 max resid 0.2302704
Run 3 stress 0.2250828
... Procrustes: rmse 0.009539203 max resid 0.1512118
Run 4 stress 0.2250839
... Procrustes: rmse 0.009650024 max resid 0.1525654
Run 5 stress 0.2250832
... Procrustes: rmse 0.009612515 max resid 0.1518641
Run 6 stress 0.2248224
... New best solution
... Procrustes: rmse 0.00631171 max resid 0.08263297
Run 7 stress 0.2253967
Run 8 stress 0.225527
Run 9 stress 0.2248128
... New best solution
... Procrustes: rmse 0.0004834727 max resid 0.004762261
... Similar to previous best
Run 10 stress 0.2255447
Run 11 stress 0.2249136
... Procrustes: rmse 0.00484646 max resid 0.08193293
Run 12 stress 0.2500137
Run 13 stress 0.2256121
Run 14 stress 0.2254928
Run 15 stress 0.2258491
Run 16 stress 0.2248191
... Procrustes: rmse 0.0008375058 max resid 0.008573354
... Similar to previous best
Run 17 stress 0.225849
Run 18 stress 0.225634
Run 19 stress 0.2249543
... Procrustes: rmse 0.006228313 max resid 0.08122835
Run 20 stress 0.2250713
... Procrustes: rmse 0.008311559 max resid 0.1485133
Run 21 stress 0.2258052
Run 22 stress 0.2248935
... Procrustes: rmse 0.004390114 max resid 0.08026974
Run 23 stress 0.2257226
Run 24 stress 0.2253296
Run 25 stress 0.2257815
Run 26 stress 0.2254095
Run 27 stress 0.2257822
Run 28 stress 0.2248496
... Procrustes: rmse 0.004509314 max resid 0.0818388
Run 29 stress 0.2250698
... Procrustes: rmse 0.008236929 max resid 0.1489688
Run 30 stress 0.225152
... Procrustes: rmse 0.01050374 max resid 0.1503719
Run 31 stress 0.2426061
Run 32 stress 0.2253598
Run 33 stress 0.418577
Run 34 stress 0.2254408
Run 35 stress 0.2254011
Run 36 stress 0.2257955
Run 37 stress 0.2256214
Run 38 stress 0.22512
... Procrustes: rmse 0.009424238 max resid 0.1483398
Run 39 stress 0.2255057
Run 40 stress 0.225153
... Procrustes: rmse 0.0106079 max resid 0.1526557
*** Best solution repeated 2 times
stressplot(df1s.mds2, main="Shepard plot")
df1s.mds2
Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE, maxit = 200)
global Multidimensional Scaling using monoMDS
Data: df1s.bc
Distance: bray
Dimensions: 2
Stress: 0.2248128
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: scores missing
3D model fits better than 2D. Stress for 2D above 0.2
Start with plots where symbol colour indicates levels of one of the factors
a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:4)],a) #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=HIV) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = HIV) )+
geom_point()+
stat_ellipse()+
labs(y="MDS3", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="HIV status") & theme_qk() & scale_color_uchicago()
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=candida) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = candida) )+
geom_point()+
stat_ellipse()+
labs(y="MDS3", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Candida") & theme_qk() & scale_color_uchicago()
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=gender) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = gender) )+
geom_point()+
stat_ellipse()+
labs(y="MDS3", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Sex") & theme_qk() & scale_color_uchicago()
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=smoking_current) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = smoking_current) )+
geom_point()+
stat_ellipse()+
labs(y="MDS3", x="MDS1")+
scale_shape_manual(values=c(16,17),
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Current smoking") & theme_qk() & scale_color_uchicago()
Try some visualizations of pairs of factors using MDS1 & MDS2 (while 3D is preferred, stress of 3D is around .15, vs 0.22, so not huge difference) Separate HIV+/-, then use symbol colour to separate second factor.
p1<-ggplot(data=subset(a, candida=="Yes"), aes(x=NMDS1, y=NMDS2, color=HIV) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title="Candida")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=subset(a, candida =="No"), aes(x=NMDS1, y=NMDS2, color = HIV) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title = "No Candida")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()
p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)
p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined &
xlim(min(p_ranges_x), max(p_ranges_x)) &
ylim(min(p_ranges_y), max(p_ranges_y))
p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color=candida) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title="HIV+")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = candida) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title = "HIV-")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()
p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)
p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined &
xlim(min(p_ranges_x), max(p_ranges_x)) &
ylim(min(p_ranges_y), max(p_ranges_y))
p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title="HIV+")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
geom_point()+
stat_ellipse()+
labs(y="MDS2", x="MDS1",
title = "HIV-")+
scale_shape_manual(
guide =
guide_legend(label.theme = element_text(size=6),
title=NULL)
)
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()
p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)
p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined &
xlim(min(p_ranges_x), max(p_ranges_x)) &
ylim(min(p_ranges_y), max(p_ranges_y))